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Abstract 

We discuss efficient conversion algorithms for orthogonal polynomials. We describe 
a known conversion algorithm from an arbitrary orthogonal basis to the monomial 
basis, and deduce a new algorithm of the same complexity for the converse operation. 
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1 Introduction 



Let (ai)j>i, (fej)j>i and (cj)j>i be sequences with entries in a field K. We can 
then define the sequence (-Fj)j>o of orthogonal polynomials in K[x] by F_i = 0, 
Fq = 1 and for z > 1 by the second order recurrence 

Fi = {a,x + hi)Fi_i + CiF,_2. (1) 

Following standard conventions, we require that OjCj is non-zero for alH > 1; 
in particular, Fi has degree i for all i > and (-Fi)j>o forms a basis of the 
K- vector space K[x]. 
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Basic algorithmic questions are then to perform efficiently the base changes 
between the basis (-Fj)i>o and the monomial basis (a;*)j>o. More precisely, for 
n E N \ {0}, we study the following problems. 

Expansion Problem (Expand^). Given ao, . . . , a^-x G IK, compute the co- 
efficients on the monomial basis of the polynomial A defined by the map 

n-l 

[ao, . . . , Qin-i] ^ A^Y^ UiFi (2) 

1=0 

Decomposition Problem (Decomp^j). Conversely, given the coefficients of 
A on the monomial basis, recover the coefficients a^, . . . , q;„_i in the decom- 
position (2) of ^4 as a linear combination of the F^'s. 

For i,j > 0, let Fij be the coefficient of a;* in Fj, and let F„ be the n x n 
matrix with entries [Fij]o<i,j<n- Problem Expand^ amounts to multiplying the 
matrix F„ by the vector [ckq, . . . , Q;„_i]*; hence, the inverse map Decomp^ is 

well-defined, since F„ is an upper-triangular matrix whose i-th diagonal entry 
Fi^i = aia2 ■ ■ ■ flj is non-zero. As we will see, the dual problem (multiplying the 
matrix F^ by a vector), denoted Expand^, plays an important role as well. 

Naive algorithms work in complexity O(n^) for both problems Expand^ and 
Decomp„. Faster algorithms are already known, see details below on prior 
work. The only new result in this article is the second part of Theorem 1 
below; it concerns fast computation of the map Decomp^. 

As usual, we denote by M a multiplication time function, such that polynomials 
of degree less than n in K[x] can be multiplied in M(n) operations in K, when 
written in the monomial basis. Besides, we impose the usual super-linearity 
conditions of [10, Chap. 8]. Using Fast Fourier Transform algorithms, M(n) 
can be taken in 0(n log(n)) over fields with suitable roots of unity, and in 
0(nlog(n) loglog(n)) over any field [18,4]. 

Theorem 1 Problem,s Expand^ and Decomp^ can be solved in 0{M(n)log{n)) 
arithmetic operations in K. 

The asymptotic estimates of Theorem 1 also hold for conversions between any 
arbitrary orthogonal bases, using the monomial basis in an intermediate step. 
In conjunction with FFT algorithms for polynomial multiplication. Theorem 1 
shows that all such base changes can be performed in nearly linear time. 

Previous work. Fast algorithms are known for problems closely related 
to Problem Expand^. From these, one could readily infer fast algorithms for 
Problem Expand^ itself. 
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In [17], the question is the computation of the values 



n— 1 

[ao, . . . , Qin-l] ^ \ Yl Oti^iiXj) ]^^.^ , 
i=0 JO<J<n 



where the the Chebyshev points xj = cos(j7r/(n — 1)). This is done by 

expanding J27=o ^i^i the Chebyshev basis and applying a discrete cosine 
transform. The article [6] studies the transposed problem: computing the map 



n-1 

[ao, . . . , an-i] ^ [ ]^^.^^- (3) 



The algorithm in [6] is (roughly, see [17] for details) the transpose of the 
one in [17]: it applies a transposed multipoint evaluation, then a transposed 
conversion, to either the monomial or the Chebyshev basis. 

Regarding problem Decomp„ to the best of our knowledge, no 0(M(n) log(n)) 
algorithm has appeared before, except for particular families of polynomials, 
like Legendre [9], Chebyshev [16] and Hermite [15]. In the case of arbitrary 
orthogonal polynomials, the best complexity result we are aware of is due 
to Heinig [13], who gives a 0(M(n) log^(n)) algorithm for solving inhomoge- 
neous linear systems with matrix F*jF„. From this, it is possible to deduce an 
algorithm of the same cost for Problem Decomp^. 



In [17], one sees mentions of left and right inverses for the related problem 

n-1 

0<j<2n-l 



n—L 



i=Q 



In [15], the inverse of the map (3) is discussed: when (xj) are the roots of F„, 
Gauss' quadrature formula shows that this map is orthogonal, so that inversion 
reduces to transposition. In other cases, approximate solutions are given. 

The various algorithms mentioned up to now have costs 0(M(n) log(n)) or 
0(M(n)log2(n)). In [2], we give algorithms of lower cost 0(M(n)) for many 
classical orthogonal polynomials (Jacobi, Hermite, Laguerre, ...), for both 
Problems Expand^ and Decomp^. 



Main ideas. Here is a brief description of the strategy used to obtain the 
complexity estimate of Theorem 1. The complete treatment with detailed 
algorithms is given in Sections 2 and 3. Three main ingredients are used: 
(i) a 0(M(n) log(n)) algorithm for Problem Expand^; (ii) the transposition 
principle; (iii) the Favard-Shohat theorem. 

We first recast (1) into the matrix recurrence [Fi,Fj_|_i]* = MW(a;)[Fj_i, Fj]*, 
where M^*) is a 2 x 2 polynomial matrix. Problem Expand^ then amounts to 
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computing T = aoM^") +aiM(^)M(°) + • • • + a„_iM("-^ • • • M(o). This is done 
by using a dividc-and-conquer algorithm similar to the one in [11, Th. 2.4] for 
the conversions between Newton and monomial bases. Assuming for simplicity 
that n is even, we rely on the decomposition T — Tq + Ti M*^5~^) • • • M(°\ with 

To = a;oM(°) + • • • + an^M^t-i) . . . m(°^ 
Ti = ctnM^t) + • • • + q;„_iM("-i) • • • M(i). 

In Section 2, a slightly different but more efficient version of this algorithm is 
given. 

An algorithmic theorem called the tmnsposifAon principle [3, Th. 13.20] states 
that the existence of an algorithm of cost 0(M(n) log(n)) for Expand^ implies 
the existence of another one with the same cost for the dual problem Expand^. 
We use an effective version of the principle, allowing to design the transposed 
algorithm in a straightforward manner starting from the direct one. 

Then, the Favard-Shohat theorem [7,19] ensures the existence of an inner 
product ( , ) on the space K[x] with respect to which the sequence is 
an orthogonal basis. This implies the matrix equality F*^ H„ F„ = D„, where 
Hn = [hi,j]o<i,j<n is the Gram matrix with hij — {x\ x^) and D„ is an invert- 
ible diagonal matrix. Its equivalent form F~^ = D~^F^H„ shows that, once 
H„ and D„ are determined. Problem Decomp„ amounts to the computation 
of the map w e K" i— > F^w G K", that is, to solving Expand^. Finally, a con- 
structive version of the Favard-Shohat theorem shows that determining the 
Gram matrix H„ can be reduced to two instances of Problem Expand^. 

In summary, by the Favard-Shohat theorem, Decomp„ is reduced to Expand„ 
and Expandjj, which can be solved in 0(M(n) log(n)), by a direct divide-and- 
conquer algorithm for the first and the transposition principle for the second. 



2 Expansion Problem 

We first describe the conversion from the orthogonal basis to the monomial 
one, and its transpose. The content of this section is mostly already known. 
However, our algorithm for the inverse operation rests crucially on these con- 
versions, so we prefer to make them explicit. 

In the following, we always suppose for simplicity that the number of unknown 
coefficients n is a power of two. For a polynomial F of degree less than m, 
m > 1, we denote by rev(F, m) = x^~-^F{l/x) the reversal of F. 
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Expansion from an orthogonal beisis. Given ao, . . . , q;„_i, we compute 
here the expansion on the monomial basis of 

A = aoFo H h an-iFn-i. 

The ideas are classical; our presentation is taken from [17]. However, our use of 
"classical" fast multiplication techniques avoids the need of precomputed con- 
stants arising in [17], and holds over any field. For i > 0, define the transition 
matrix 

^ 1 



so that we have 



For J > i, let M^^J) = M^^-^'^^M^^-^'J-i) • . . m(^'^+^); for i = j, M^^J) is the 
2x2 identity matrix. It follows that we have 



Fi 


— ]y[(*,*+i) 











F,-i 




F.-i 




_F, _ 







besides, for i > j > i, we have the associativity relation M*^*'' 
We can then rewrite A as 



Fo 
Fi 



H h 



an-2 OLn-\ 



Fn-2 
Fn-1 



where the sum has n/2 terms. We deduce the equalities 



A 



Fo 
Fi 



+ ••• + 



Oin-2 CXfi-l 



Fo 




Fo 


= B 


Fi_ 




Fi 



n/2-l 



where B is the 1x2 matrix B = 



Oi2i Oi2i+l 



The computation of A is thus reduced to that of the matrix B. Write n' — n/2. 
Following [20] and [14], we build the subproduct tree associated to the transition 
matrices M'^-''^^ This is a complete binary tree having d — log2(n) — log2(n')-|-l 
rows of nodes labeled as follows: 

• the leaves of the tree are labeled by the matrices L*^*^"^'*) = ]V[(2i+i,2i+3)^ f^j. 
i — 0, . . . ,n' — 1; 
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Expand(a, b, c, A, n) 



Input: A = Y^i^Q QiX^ and a, b, c 
Output: Er=o OiiFi 



L(:/» ^ SubproductTree(a, b, c, n) 
for i = 0, . . . , 2*^-1 - 1 do 



(^1 ^ «2j+l 

for j = - 2, . . . , do 




for z = 0, . . . , 2^ - 1 do 




(i+l,2i+l) ;-(i+l,2j) 

-^0,0 

(i+l,2i+l) j{j+l,2i) 

-^0,1 



+ y 



(i + l,2i + l) ;-(j + l,2i) 

1 -^1,0 

(j + l,2i + l) r(j + l,2i) 

1 -^1,1 



Fig. 1. Algorithm solving Problem Expand, 



n 



• for J = 0, . . . , (i — 2, there are 2^ nodes of depth j and the (1 + i)-th one is 
indexed by the matrix L^^'^) = L(j+i,2i+i)L(j+i,2i)^ for < i < 2^' - 1. 

The entries /^^-'^J-' of L*^-''*) have degrees at most 2'^~^ with < m, f < 1. 
An easy induction also shows that for j = 0, . . . , d — 1 and i = 0, . . . , 2-^ — 1, 
we have the equality 



The cost of computing all matrices in the tree is 0(M(n) log(n)), as in [10, 
Chapter 10]. Then, to compute B, we go up the subproduct tree and perform 
linear combinations along the way: we maintain a family of 1 x 2 vectors 



^ ^p-.O]^ ^j^jj j = 0, . . . , d - 1 and i = 0, . . . , 2^' - 1, such that 



The overall cost is again 0(M(n) log(n)). 

Remark that not all the nodes of the complete subproduct tree are actually 
needed in this algorithm. Indeed, its rightmost branch containing L*^-''^^"^^ for 
< i < <i — 1 is not necessary in the computation described in Equation (4). 

In the pseudo-code in Figure 1, we make all scalar operations explicit, so 

as to make the transposition process easier in the next paragraph. Starting 
from the sequences a = (ai, . . . , a„_i), b = . . . , 6„_i), c = (ci, . . . , c„_i), 
the subroutine Subproduct Tree(a, b, c, n) computes the matrices M^*'*"*"^) for 
< i < n - 2, then the matrices L^^'*) for 1 < j < d - 1 and < i < 2^' - 2. 
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Expand*(a, b, c, A, n) 

Input: A = Y^^Iq Q.iX'^ and a, b, c 
Output: ( YJlZo a^coeff (Fj-, i) )j=o,...,n-i 

LO.*) ^ SubproductTree(a, b, c, n) 



(0,0) 



.(0,0) 



mul*(/l,Fo,5rf,o) 
mul*(A,Fi,5;,,o) 



for j = 0, . . . , - 

for i = 2^ - 1, 

.,(i+i,2i) 



2 do 

. . , do 

-Vq'^^ mod 



(j+l,2i+l) 
^0 



rj+i,2i+i) 



mod 

^mul*(4^'^L^o'-'''^<^.,.+i) 
+ mul*(^P\LS:,+^'^^<5,,+0 

+ mul*(^i^'^Li^^'^^5',^,) 



return 



(d-1,0) (rf-1,0) 



d-l 



-1) (d-l,2'*-i-l) 



Fig. 2. Algorithm solving Problem Expand^ 

Transposed expansion. Let r, s > 1 and let M be a r x s matrix with 
entries in K. The transposition principle [3, Th. 13.20] states that the existence 
of an algorithm for the matrix-vector product h ^ M6 implies the existence 
of an algorithm with the same cost, up to 0(r + s) operations, to perform 
the transposed matrix-vector product c h- > M*c. This paragraph gives the 
transposed version of the conversion algorithm above: a similar algorithm is 
given in [6], but our derivation is substantially more compact. 

A fundamental operation is transposed polynomial multiplication. For k in N, 
let 'K[x\k be the K- vector space of polynomials of degree less than k. Then, for 
B in ]K[a;] of degree m, we let mul(., k) be the multiplication-by-S operator, 
defined over Kfx]^; its image lies in 

The transpose of this map is denoted by mul*(., B, k)] by identifying ^[x\k with 
its dual, one sees that mul*(.,S,A;) maps Kfxjfc+m to K[x]fe. In [1,12], details 
of the transposed versions of plain, Karatsuba and FFT multiplications are 
given, with a cost matching that of the direct product. Without using such 
techniques, writing down the multiplication matrix shows that mul*(., B,k) is 

A e K[x\k+m ^ {A rev(B, m + 1) mod div x"" e K[x]k. 

Using standard multiplication algorithms, this formulation leads to slower 
algorithms than those of [1,12]. However, here k and m are of the same order 
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of magnitude, and only a constant factor is lost. 

Using this tool, the transposed expansion algorithm in Figure 2 is obtained 
by reversing the flow of the direct one in Figure 1. The loops are traversed 
in opposite order. Then, the operation v^-''*) = v(-^+^'^*) + v(-^+^'^*+^)L(-'+^'^*^ 
in the inner loop is replaced by a truncated copy of v*^-''*^ into v*^-'+^'^*) and 
a transposed matrix-vector product, where polynomial multiplications arc re- 
placed by transposed multiplications. To perform truncations and transposed 
multiplications, we need information on the degrees of the polynomials in- 
volved. By induction, we get the following inequalities, ior j — 0, . . . ,d — 1 
and i = 0, ...,2^' - 1, 

deg(4^''^) < Sdj = max(0, 2<^-^' - 3) + 1, deg(t;p'^^) < 6'^^^ = 2'^-^ - 1. 

This information enables us to write the transposed algorithm in Figure 2. 
Using either the transposition principle or a direct analysis, one sees that the 
cost of this algorithm is 0(M(n) log(n)). 



3 Decomposition Problem 

The Favard-Shohat theorem [7,19], see also [5, Theorem 4.4], asserts that for 
(Fj) as in (1), there exists a linear form L : K[x] K for which (Fj) is formally 
orthogonal, in the sense that, for i > 1, 

L{FiFj) = Q for < j < i, L{Fl)^Q. 

The linear form L is specified by its moments L{x^), for i > 0, or equivalently 
by the generating series 

For completeness, we give in the following theorem a self-contained, construc- 
tive, proof of this classical result, showing how to compute truncations of Sl- 
The proof is inspired by the presentation in [8, Section 3]. 

Theorem 2 Let (Fi) be the sequence satisfying F^i = 0, Fq = 1 and recur- 
rence (1). Define the sequence (Gi) by G-i = 0, Gq = 1 and, for i > 1 

Gi — {tti-^ix + bi+i)Gi-i + Ci+iGi-2- 
Then, there exists a K-linear form L : K[x] — > K such that 
LiFiFj)^0 for iy^j, and L{F^) ^ {-ly " " " ^'+^ for i>0. (5) 
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Moreover, for any i > 1, the following equality holds between truncated series 
in K[[x]]: 

revlG,-!,^^^^^^,^^, modx2\ (6) 



rev(Fi,i + l) 



i>0 



Proof. For i > 0, write F* = rev(Fj, i + 1) and G* = rev{Gi, i + Let also 
define F*^ — G*_-^ — 0. These polynomials satisfy the recurrences 

F* = {ai + bix)F,U + CiX^F*_^, G* = (a^+i + bi+ix)Gti + Ci+ix^G*_2, 
for i > 1, which can be recast into the matrix form 



f: gu 









1 



F* GU 



Taking determinants, we deduce that for ^ > 1 the following identity holds 



Fi+i 



= -Ci+l 



TP* / r~1-k 



G* 



F.^ 



-X 



i-2 



i+1 



FU 



Applying it to i, i — 1, . . . and denoting 7, = C2 • • • Q, we get that for i > 1, 

G1 



F*+i 



F* 



z?* z?* 



(7) 



A separate check shows that Equation (7) also holds for i = 0. 



For i > 0, F* has constant coefficient 5i — Oi • • • Oj, which is non-zero, and 
is thus invertiblc in IK[[a;]]. Since the 7^+1 are non-zero as well. Equation (7) 
shows that the sequence G*/F*^i is Cauchy and thus convergent in IK[[a;]]. 
Besides, if we let S be its limit, summing up Equation (7) for i, i + 1, . . . yields 



+ i-iy^^x''' modx^'+\ for i>0. 

OiOi+i 



(8) 



Write S = J: i>o ^iX^ and define the linear form L on K[x] by L(x*) = £j. Then 
Equation (6) is a direct consequence of (8). 

For i > 0, equating coefficients of x*, . . . , x^*""^ and in Equation (8) multi- 
phed by F* implies L{FiX^) = for i < j and L{FiX') = (-l)'7j+i/(5i+i. By 
hnearity, this shows that L also satisfies Equality (5). □ 



Proof of Theorem 1. We can now prove the second part of Theorem 1, 
dealing with expansions in the monomial basis. The corresponding algorithm 
is given in Figure 3. 
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Decomp(a, b, c, A, n) 

Input: A ~ YhZq UiX^ and a, b, c 

Output: ao,..., q;„_i such that A — J2i=o <^iFi 

a! cat(a, 1) 
b' ^ cat(b, 1) 
c' cat(c, 1) 

F <— Expand(a', b, c, x", n+1) 

G ^ Expand(Sa', Sb', Sc', n) 

Q ^ rey{G, n) /rev(F, n + 1) mod x^^-^ 

y ^mul*(g,A n) 

w Expand^ (a, b, c, V, n) 

di <— (-1)^2 • • -Ci+i/ai+i, ioT <i <n 

ai <— Wi/di for < i < n, where w — [wq, . . . , Wn-iY 

return ckq, . . . , an-i 

Fig. 3. Algorithm solving Problem Decomp„ 

We first compute (L(a;*))j<2n-i- To do this, we start from the sequences a, b 
and c to which we add the element 1, in order to make the polynomial F = 
well-defined (any non-zero choice would do). We then use the algorithm Ex- 
pand of the previous section to compute G = G^-i and F — F^ and we 
determine the power series expansion rev(G'„_i, n)/rev(F„, n + 1) mod x^""^. 
The first step takes 0(M(n) log(n)) operations, and the second one 0(M(n)) 
using Newton iteration [10, Chap. 9]. In the pseudo-code we use the nota- 
tion Sx for the shifted sequence (xj+i) of x = (xj) and the notation cat for 
concatenation. 

Consider finally the matrix F„ defined in the introduction, and let H„ = 
[ifij]o<i,j<n be the nxn Hankel matrix with iJ^j = L{x''^^). Let next D„ be 
the diagonal matrix of size n, with — L{Ff). We deduce the factorization 

F^H„F, = D,, or F;^ = D;;^ F^ H„. 

Equation (5) shows that one can compute the entries of D„ in 0{n) operations. 

At this stage, all elements of D„ and H„ are known. Right-multiplication 
of H„ by the coefficient vector of a polynomial A e ]K[a;]„ amounts to the 
transposed multiplication of the polynomial Q — Z^^"o ^ L{x^)x^ by A, that can 
be performed in time M{n) +0{n). Using the transposed expansion algorithm 
Expand* of the previous section, multiphcation by F^ costs 0(M(n) log(n)). 
Finally, multiplying by takes linear time. This concludes the proof of 
Theorem 1. 
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